package scu.maqiang.homogenization;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L2PE;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

/**
 * 平面轴对称压电特征值的二阶双尺度有限元计算
 * @author 马强
 *
 */
public class AxisPEEigenValueSOTS1D {

	public static void main(String[] args) {
		
		String dir = "AxisPEEigenValue\\";
		int N = 32;
		double ep = 1.0 / N;

		ScalarFunc coefG1 = (xyz, label, param) -> label == 10 ? 200 : 4;
		ScalarFunc coefG2 = (xyz, label, param) -> label == 10 ? 300 : 6;
		ScalarFunc coefLambda = (xyz, label, param) -> label == 10 ? 150 : 3;
		ScalarFunc coefEp1 = (xyz, label, param) -> label == 10 ? 1.0 : 0.02;
		ScalarFunc coefE1 = (xyz, label, param) -> label == 10 ? 5.0 : 0.3;
		ScalarFunc coefE2 = (xyz, label, param) -> label == 10 ? 6.0 : 0.4;
		ScalarFunc coefrho = (xyz, label, param) -> label == 10 ? 1 : 0.05;
		ScalarFunc coeff1 = (xyz, label, param) -> label == 10 ? 1 : 1;
		ScalarFunc coefrf = (xyz, label, param) -> label == 10 ? 1 : 1;

//		Mesh1L2Generator mg = new Mesh1L2Generator();
		Mesh1L2 cellMesh = new Mesh1L2();
		cellMesh.line(0, 1.0, 256);
		ScalarFunc cellRegionFunc1D = (x, llable, param) -> Math.abs(x[0] - 0.5) < 0.25 ? 1.0 : -1.0;
		cellMesh.setDomainLabel(cellRegionFunc1D, null, 10);

		FES1L2PE fsCell = new FES1L2PE(cellMesh);
		double[][] FOCS = new double[3][fsCell.GetNdof()];
		double[][] SOCS = new double[6][fsCell.GetNdof()];
		double[][] D = new double[3][3];
		SOTS1DL2PE ss = new SOTS1DL2PE(fsCell);
		ss.AxisPECellSolution(new ScalarFunc[] { coefG1, coefG2, coefLambda, coefEp1, coefE1, coefE2, coefrho, coeff1, coefrf },
				null, FOCS, D, SOCS);
		double[][] FOCSUV1 = new double[2][cellMesh.getNv()];
		double[][] FOCSUV2 = new double[2][cellMesh.getNv()];
		double[][] FOCSUV3 = new double[2][cellMesh.getNv()];
		fsCell.extractComponent(FOCS[0], FOCSUV1);
		fsCell.extractComponent(FOCS[1], FOCSUV2);
		fsCell.extractComponent(FOCS[2], FOCSUV3);
		double[][] FFOCS = new double[][] { FOCSUV1[0], FOCSUV1[1], FOCSUV2[0], FOCSUV2[1], FOCSUV3[0], FOCSUV3[1] };
		System.out.println("N1 min: " + MVO.min(FOCSUV1[0]) + "\tmax: " + MVO.max(FOCSUV1[0]));
		System.out.println("H  min: " + MVO.min(FOCSUV1[1]) + "\tmax: " + MVO.max(FOCSUV1[1]));
		System.out.println("F  min: " + MVO.min(FOCSUV2[0]) + "\tmax: " + MVO.max(FOCSUV2[0]));
		System.out.println("R1 min: " + MVO.min(FOCSUV2[1]) + "\tmax: " + MVO.max(FOCSUV2[1]));
		System.out.println("M  min: " + MVO.min(FOCSUV3[0]) + "\tmax: " + MVO.max(FOCSUV3[0]));
		System.out.println("S  min: " + MVO.min(FOCSUV3[1]) + "\tmax: " + MVO.max(FOCSUV3[1]));

		cellMesh.toTecplot(dir + "CellFOCS.dat", FFOCS);

		Mesh1L2 homoMesh = new Mesh1L2();
		homoMesh.line(0, 1, 2048);
		homoMesh.translate(1.0);
		FES1L2PE fsHomo = new FES1L2PE(homoMesh);
		SRMatrix homoA = new SRMatrix(fsHomo.GetNdof());
		SRMatrix homoM = new SRMatrix(fsHomo.GetNdof());
		double[] homoX = new double[fsHomo.GetNdof()];
		fsHomo.assemblePEStiff(new double[] { D[0][0], D[0][1], D[0][2], D[1][0], D[1][1], D[1][2] },
				BVPType.AXISYMMETRIC, homoA);
		fsHomo.assemblePEMass(new double[] { D[2][0] }, BVPType.AXIS_CONSIST_MASS, homoM);
		Arrays.fill(homoX, 1.0);
		// fsHomo.applyBC_MBN(homoA, homoX, Direct.X, 0.0, 1);
		// fsHomo.applyBC_MBN(homoA, homoX, Direct.Y, 0.0, 2);

		int numEigen = 10;
		double[] homoD = new double[numEigen];
		double[][] homoV = new double[numEigen][fsHomo.GetNdof()];
		EigenSSolver eigenSolver = new EigenSSolver();
		eigenSolver.inversePowerTran_PLU(homoA, homoM, fsHomo.getnBand(), homoX, 20.0, numEigen, homoD, homoV);

		Mesh1L2 fineMesh = new Mesh1L2();
		fineMesh.line(0, 1, 2048);
		ScalarFunc regionFineFunc1D = (x, llable, param) -> {
			double xCell = x[0] * N - (int) (x[0] * N);
			return Math.abs(xCell - 0.5) < 0.25 ? 1.0 : -1.0;
		};
		fineMesh.setDomainLabel(regionFineFunc1D, null, 10);
		fineMesh.translate(1.0);

		FES1L2PE fsFine = new FES1L2PE(fineMesh);
		SRMatrix fineA = new SRMatrix(fsFine.GetNdof());
		SRMatrix fineM = new SRMatrix(fsFine.GetNdof());
		double[] fineX = new double[fsFine.GetNdof()];
		fsFine.assemblePEStiff(new ScalarFunc[] { coefG1, coefG2, coefLambda, coefEp1, coefE1, coefE2 }, null,
				BVPType.AXISYMMETRIC, fineA);
		fsFine.assemblePEMass(new ScalarFunc[] { coefrho }, null, BVPType.AXIS_CONSIST_MASS, fineM);
		Arrays.fill(fineX, 1.0);
		// fsFine.applyBC_MBN(fineA, fineX, Direct.X, 0.0, 1);
		// fsFine.applyBC_MBN(fineA, fineX, Direct.Y, 0.0, 2);

		double[] fineD = new double[numEigen];
		double[][] fineV = new double[numEigen][fsFine.GetNdof()];
		eigenSolver.inversePowerTran_PLU(fineA, fineM, fsFine.getnBand(), fineX, 20.0, numEigen, fineD, fineV);

		double[] D1 = new double[numEigen];
		double[] D2 = new double[numEigen];
		double[][] V0 = new double[numEigen][fsFine.GetNdof()];
		double[][] V1 = new double[numEigen][fsFine.GetNdof()];
		double[][] V2 = new double[numEigen][fsFine.GetNdof()];
		// ss.assemblePEEigenSOTSSolution(fsFine, V0, V1, V2, D1, D2, fsHomo, homoV,
		// homoD, fsCell, FOCS, SOCS, coefrho, N, ep);

		for (int i = 0; i < numEigen; i++) {
			ss.assembleAxisPESOTSSolution(fsFine, V0[i], V1[i], V2[i], fsHomo, homoV[i], fsCell, FOCS, SOCS, N, ep);
		}

		ss.assembleAxisPEEigenSOTSSolution(fsFine, V0, V1, V2, coefrho, homoD, D1, D2, ep);

		System.out.println("Fine      Homo    FOTS      SOTS      Error0      Error1      Error2");
		for (int i = 0; i < numEigen; i++) {
			System.out.println(fineD[i] + "\t" + homoD[i] + "\t" + D1[i] + "\t" + D2[i] + "\t"
					+ Math.abs(homoD[i] - fineD[i]) / fineD[i] + "\t" + Math.abs(D1[i] - fineD[i]) / fineD[i] + "\t"
					+ Math.abs(D2[i] - fineD[i]) / fineD[i]);
		}
		
		String outputFilePreFix = dir + "AxisPEEigenValueSOTS";
		double[][] fineU = new double[2][fineMesh.getNv()];
		double[][] U0 = new double[2][fineMesh.getNv()];
		double[][] U1 = new double[2][fineMesh.getNv()];
		double[][] U2 = new double[2][fineMesh.getNv()];
		double[][] uv012 = new double[4][fineMesh.getNv()];
		double[][] phiv012 = new double[4][fineMesh.getNv()];
		int fineNv = fineMesh.getNv();
		for (int i = 0; i < numEigen; i++) {
			fsFine.extractComponent(fineV[i], fineU);
			fsFine.extractComponent(V0[i], U0);
			fsFine.extractComponent(V1[i], U1);
			fsFine.extractComponent(V2[i], U2);

			double maxU = MVO.maxA(fineU[0]);
			double maxU0 = MVO.maxA(U0[0]);
			double maxU1 = MVO.maxA(U1[0]);
			double maxU2 = MVO.maxA(U2[0]);
			for (int j = 0; j < fineNv; j++) {
				uv012[0][j] = fineU[0][j];
				uv012[1][j] = U0[0][j] / maxU0 * maxU;
				uv012[2][j] = U1[0][j] / maxU1 * maxU;
				uv012[3][j] = U2[0][j] / maxU2 * maxU;
			}

			double maxPhi = MVO.maxA(fineU[1]);
			double maxPhi0 = MVO.maxA(U0[1]);
			double maxPhi1 = MVO.maxA(U1[1]);
			double maxPhi2 = MVO.maxA(U2[1]);
			for (int j = 0; j < fineNv; j++) {
				phiv012[0][j] = fineU[1][j];
				phiv012[1][j] = U0[1][j] / maxPhi0 * maxPhi;
				phiv012[2][j] = U1[1][j] / maxPhi1 * maxPhi;
				phiv012[3][j] = U2[1][j] / maxPhi2 * maxPhi;
			}
			fineMesh.toTecplot(outputFilePreFix + (i + 1) + "u.dat", uv012);
			fineMesh.toTecplot(outputFilePreFix + (i + 1) + "phi.dat", phiv012);
		}
	}
}
